function dz = diff_Reg3(t,z)

% Differential equations in Regime 3, where j>0, iF>=0, renewable only.

global par

k = z(1);
lambda = z(2);
kB = z(3);
h = z(4);
f = z(5);
phi = z(6);
eta = z(7);

dz = zeros(7,1);
if eta == 0
    jnv = 0;
    bhdotkB = 0;
else
    jkBratio = (lambda/eta/(1-par.psi))^(-1/par.psi);
    jnv = jkBratio * kB;
    bhdotkB = par.b * kB^(1+ par.psi) * jnv^(1-par.psi);
end
c = lambda^(-1/par.gamma);

if lambda == phi

    pe = lambda * (par.A - par.delta) / (par.a*k + par.F0 - par.a * f);
    LHS = zeros(3,3);
    RHS = zeros(3,1);

    LHS(1,:) = ones(1,3);
    LHS(2,:) = [(par.F0 - par.a*f), -(par.H0 + par.b*h), -par.a*k];
    LHS(3,1) = (par.m+par.A-par.psi*jkBratio/(1-par.psi))*a;
    LHS(3,3) = (par.m+par.delta-par.psi*jkBratio/(1-par.psi))*a;

    RHS(1) = par.A*k - par.m*kB - jnv - c;
    RHS(2) = bhotkB;
    RHS(3) = LHS(3,1)*par.delta*k + (par.A-par.delta)*par.b* ...
        jkBratio^(1-par.psi)*kB + (par.F0 + par.a*k + par.a*f) ...
        * (1/eta - 1/lambda) * jkBratio/(1-par.psi);

    sol = RHS\LHS;

    inv = sol(1);
    iB = sol(2);
    iF = sol(3);

    dz(2) = (par.beta + par.delta - par.A) * lambda + (par.F0-par.a*f) * pe;
else

    inv = ((par.H0+par.b*h) * (par.A*k - par.m*kB - jnv - c) + ...
        bhdotkB) / (par.Fbar + par.H0 + par.b*h);

    iB = (par.Fbar*inv - bhdotkB)/(par.H0+ par.b*h);

    iF = 0;

    pe = ((par.m+par.A)*lambda - par.psi*eta*(jnv/kB)^(1-par.psi)) ...
        / (par.Fbar+par.H0+par.b*h);

    dz(2) = (par.beta + par.delta - par.A) * lambda + par.Fbar * pe;
end

dz(1) = inv - par.delta*k;

dz(3) = iB - par.delta * kB;
dz(4) = kB^par.psi * jnv^(1 - par.psi);
dz(5) = iF;
dz(6) = par.beta * phi - par.a * pe * k;
dz(7) = par.beta * eta - par.b * pe * k;
